
library(sf)
library(tidyverse)
library(kableExtra)

# Loading the Shape files

brazil <- st_read("./_3_data/shapefile_municipality/gadm40_BRA_2.shp")

pede <- brazil %>% filter(NAME_2 %in% c("Pederneiras"))
soro <- brazil %>% filter(NAME_2 %in% c("Sorocaba"))

#Loading urban area shape file

urban_pede <- st_read("./_3_data/shapefile_urban/pede_urban.shp")
urban_soro <- st_read("./_3_data/shapefile_urban/soro_urban.shp")

#Creating boundary of urban area as separate object
urban_bound_pede <- st_boundary(urban_pede)
urban_bound_soro <- st_boundary(urban_soro)

# Identifying urban area within municipality:

sf::sf_use_s2(FALSE) #Required for using spherical geometry

#Intersection between urban and municipality boundaries

inter_pede <- st_intersection(st_make_valid(urban_pede),pede)
inter_soro <- st_intersection(st_make_valid(urban_soro),soro)


#Clinic coordinates

clinic_pede <- read.csv("./_3_data/ubs.csv",fileEncoding="UTF-8-BOM") %>% filter(municipality=="Pederneiras")
clinic_soro <- read.csv("./_3_data/ubs.csv",fileEncoding="UTF-8-BOM") %>% filter(municipality=="Sorocaba")

#Creating sf objective to enable calculation of distance
#Converting coordinate system to enable distance calculation

clinic_coords_pede <- st_as_sf(clinic_pede,coords=c("longitude", "latitude")) %>% st_set_crs(4326) %>% st_transform(crs=32617)
clinic_coords_soro <- st_as_sf(clinic_soro,coords=c("longitude", "latitude")) %>% st_set_crs(4326) %>% st_transform(crs=32617)

#Distance 

pede_mean <- NA
soro_mean <- NA
pede_sd <- NA
soro_sd <- NA

#Overlaying grid on the intersection map and caclulating distance

grid_spacing <- c(.001,.005,.01)

#Creating overlap 
for(i in 1:3){
overlap_pede <- inter_pede %>%
  st_make_grid(cellsize =c(grid_spacing[i],grid_spacing[i])) %>%
  st_intersection(inter_pede)  %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>% 
  st_transform(32617)

overlap_soro <- inter_soro %>%
  st_make_grid(cellsize =c(grid_spacing[i],grid_spacing[i])) %>%
  st_intersection(inter_soro)  %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>% 
  st_transform(32617)

# Centroids of grid blocks

block_centroid_pede <- overlap_pede %>%  st_centroid()
block_centroid_soro <- overlap_soro %>%  st_centroid()

#Distance of centroids to each clinic 

distance_pede <- as.data.frame(st_distance(block_centroid_pede,clinic_coords_pede))
distance_soro <- as.data.frame(st_distance(block_centroid_soro,clinic_coords_soro))

#Distance of centroids to the closest clinic 
#Divide by 1000 as output is in metres.

minimum_distance_pede <- (apply(distance_pede, 1, FUN = min))/1000 
minimum_distance_soro <- (apply(distance_soro, 1, FUN = min))/1000

#Final output:

pede_mean[i] <- mean(minimum_distance_pede)
soro_mean[i] <- mean(minimum_distance_soro)

pede_sd[i] <- sd(minimum_distance_pede)
soro_sd[i] <- sd(minimum_distance_soro)
}

Municipality <- c("Pederneiras","Pederneiras","Pederneiras","Sorocaba","Sorocaba","Sorocaba")
Degree <- c(.001,.005,.01,.001,.005,.01)
mean <- c(pede_mean,soro_mean)
sd <- c(pede_sd,soro_sd)

table <- data.frame(Municipality,Degree,distance=mean,SD=sd)

kbl(table,format = "latex")

write.table(table, file = "./_4_outputs/table_oa_iv.htm", sep=" ", row.names=TRUE, col.names=TRUE)